
% ls_scale
% Christopher J. Whalen
% University of Illinois
% 08/25/05 - Date Created
% 08/25/05 - Last Modified

% Purpose - to scale the data by fitting it to a least squares regression.
% The program will scale from the centroid

function [xyz,err] = ls_scale(xyz,xyz_mask,mri_dim);


fprintf('\nScaling Data...\n');
% stop_scale = 0.01; % Scaling stops when the largest scaling less than 1%
% j= 0; % counter for scaling iterations
% while stop_val>=
%     j = j + 1;
%     fprintf('Iteration %s\n',j);

xyz_dist = zeros(size(xyz_mask,1),1); % preallocating 
xyz_scalp_pt = zeros(size(xyz,1),3);  % preallocating
min_dist_vec = zeros(size(xyz,1),1);

% Finding Closest Point
for i=1:size(xyz,1)
    xyz_repmat = repmat(xyz(i,:),size(xyz_mask,1),1);
    xyz_diff = xyz_mask - xyz_repmat;
    xyz_dist(:) = xyz_diff(:,1).^2 + xyz_diff(:,2).^2 + xyz_diff(:,3).^2;
    [min_dist,index] = min(xyz_dist(:)); % finding the nearest point
    min_dist_vec(i,1) = min_dist;
    index = index;
    xyz_scalp_pt(i,:)=xyz_mask(index,:);
end

err = xyz_scalp_pt - xyz;

Sy1 = std(xyz_scalp_pt(:,1));
Sx1 = std(xyz(:,1));
R1 = 0;
for i=1:length(xyz)
    R1 = R1 + (xyz(i,1) - mean(xyz(:,1)))*(xyz_scalp_pt(i,1) - mean(xyz_scalp_pt(:,1)))/(Sx1*Sy1);
end
R1 = R1/(length(xyz)-1);
%R1 = corrcoef(xyz(:,1),xyz_scalp_pt(:,1))
a1 = R1*Sy1/Sx1;
b1 = mean(xyz_scalp_pt(:,1)) - a1*mean(xyz(:,1));


Sy2 = std(xyz_scalp_pt(:,2));
Sx2 = std(xyz(:,2));
R2 = 0;
for i=1:length(xyz)
    R2 = R2 + (xyz(i,2) - mean(xyz(:,2)))*(xyz_scalp_pt(i,2) - mean(xyz_scalp_pt(:,2)))/(Sx2*Sy2);
end
R2 = R2/(length(xyz)-1);
%R2 = corrcoef(xyz(:,2),xyz_scalp_pt(:,2))
a2 = R2*Sy2/Sx2;
b2 = mean(xyz_scalp_pt(:,2)) - a2*mean(xyz(:,2));


Sy3 = std(xyz_scalp_pt(:,3));
Sx3 = std(xyz(:,3));
R3 = 0;
for i=1:length(xyz)
    R3 = R3 + (xyz(i,3) - mean(xyz(:,3)))*(xyz_scalp_pt(i,3) - mean(xyz_scalp_pt(:,3)))/(Sx3*Sy3);
end
R3 = R3/(length(xyz)-1);
%R3 = corrcoef(xyz(:,3),xyz_scalp_pt(:,3))
a3 = R3*Sy3/Sx3;
b3 = mean(xyz_scalp_pt(:,3)) - a3*mean(xyz(:,3));



% hold on;
% x1 = -300:1:300;
% y1 = a1*x1 + b1;
% plot(x1,y1)
% plot(xyz_scalp_pt(:,1),xyz(:,1));
% pause;
% close;
% hold on;
% x2 = -300:1:300;
% y2 = a2*x2 + b2;
% plot(x2,y2)
% plot(xyz_scalp_pt(:,2),xyz(:,2));
% pause;
% close;
% hold on;
% x3 = -300:1:300;
% y3 = a3*x3 + b3;
% plot(x3,y3)
% plot(xyz_scalp_pt(:,3),xyz(:,3));
% pause;
% close;


% apply scaling
for i=1:length(xyz)
    xyz(i,1) = a1*xyz(i,1) + b1;
    xyz(i,2) = a2*xyz(i,2) + b2;
    xyz(i,3) = a3*xyz(i,3) + b3;
end

